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FORTRAN PROGRAM FOR COMPUTATION OF WEBER 
FUNCTIONS AND THEIR FIRST DERIVATIVES 


by Antra Priede and Gabriel Allen 
Lewis Research Center 

SUMMARY 

This report contains a description of a FORTRAN IV program for computing D^x) 
and (d/dx)D j/ (x), which are the Weber function of order v and argument x and its first 
derivative, respectively. The range of values is x[0, 30] and i[ 1,200]. The computa- 
tions are performed in double precision and the values are accurate to seven significant 
figures over most of this region. However, in the range x[5.5, 50] and (x 2 /4) + 1 

< v < (x 2 /4) +13, the accuracy is reduced. Included herein are tables of D^(x) and 
(d/dx)D l/ (x) with values listed at intervals of 1 in x and 5 in v; the tables begin with 
x = 0 and v = 4. 5, but they omit part of the region covered by x[5. 5, \AjO]. 

The program is so written to take advantage of asymptotic forms, recursion relations^ 
and series computations at whatever combinations of v and x these methods yield the 
desired accuracy in the shortest machine time. 


INTRODUCTION 

Many problems of physical and engineering interest are concerned with waves ema- 
nating into and out of systems having parabolic geometry (ref. 1). Separation of the ap- 
propriate wave equation (either classical or quantum -mechanical) in parabolic coordinates 
results in a differential equation that can be transformed into Weber’s equation. The 
latter may be written as 


fv + I 
dx 2 V 2 


H’- 


( 1 ) 


The solutions to this equation are called Weber functions or parabolic cylinder func- 
tions. They are commonly denoted by D^(x), where both x and v may take on general 



complex values. Physically, x usually signifies a distance from the axis of the parabolic 
region, while the order v contains other physical properties intrinsic to the system. 
Thus, the solutions of greatest interest are those for which both x and v are real. 

For physical applications it is convenient to use, as independent solutions to Weber's 
equation (1), functions distinguished by their asymptotic behavior. Thus, D y (x) goes to 0 
as x — whereas a second solution, which is linearly independent of D [/ (x), approaches 
infinity as x — °°. There is an important connection between this second solution and the 
form of the first solution (i. e. , of D^(x)) in the range x < 0. 

As indicated in reference 2, Weber's equation is defined for complex arguments z. 
The asymptotic expansion of D v (z) is valid in the sector of the complex plane given by 
|arg z | < 37r/4. Note that the negative part of the real axis (arg z = n) is not in this 
region of validity. Thus, the use of the solution to Weber's equation which is valid for 
x < 0 will require the second independent solution. This solution may take any of the 
forms D_j,_j(i z), D_ l ,_j(-iz), or D J/ (-z). These three forms are not identical, but there 
are two features common to each of them: first, the region of validity includes 
arg z = 7r; and second, the forms are linearly independent of D^(z). 

For 5 tt/ 4 > arg z > 7t/4, 

D „ .(-te) 

r(-i') 


and for -7 t/ 4 > arg z > - (57 t/4), 

D (z) ~ e~ Uni D (-z) + j/ll"e"^ +1 ^ 7ri//2 D 1 (iz) 

v v r(-v) ~ v ~ x 

where the symbol ~ means approaches asymptotically. 

Because of these properties, one should not substitute negative values of the argu- 
ment into D^(x) before carefully examining reference 2 in order to obtain the proper 
combination of second solutions which have the behavior desired. 

This particular program was developed in order to solve a quantum -mechanical 
problem in solid-state physics (ref. 3) by using solutions of Weber's equation for a dense 
distribution of v in the interval [1, 200] at a fixed value of x (about 12.5). The most 
complete tables and discussion of D^(x) which are currently available (ref. 4) cover the 
range v[-5. 5, 4. 5] and x[0, 5]. The use of these tables for the range of orders required 
would necessitate repeated recursion (see eq. (14)) followed by interpolation. Further- 
more, the desired value of x lay outside the covered range. For these reasons, it was 
necessary to devise a method for computing D^x) for the desired values of v and x. 

A FORTRAN IV program which could accomplish this task was written for the Lewis 
7094II - 7044 direct-couple system. This report describes a modified program which 
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can be used to evaluate D^(x) and (d/dx)D^(x) in the range x[0, 30] and v[ 1, 200]. Com- 
putations are made only for that solution to equation (1) which goes to zero as x — 


ANALYSIS OF COMPUTATIONAL PROCEDURES FOR THE WEBER FUNCTION 


In what follows, the notation of reference 4 will be used. The relations between this 
notation and the definitions in the introduction are 

a = -v - — (2) 

2 

U(a,x) = D^x) (3a) 

In this connection, the following relation will be used in Order to compute (d/dx)U(a,x): 

— U(a, x) = - - U(a, x) - (a + -Vj(a+1, x) (3b) 

dx 2 \ 2/ 

Weber's equation in this notation becomes 



Reasonable procedures for computing values of U(a,x) and its derivative with an 
accuracy of seven significant figures are presented herein. The range covered is 
x[0, 30] and a[-l. 5, -200. 5]. Thus, the properties described in this section are only 
those which are relevant to this objective. As might be expected, over the range of a 
and x desired, several different methods of computation must be used to obtain the re- 
quired accuracy efficiently. Discussion of this matter will be facilitated by reference to 
figure 1 where that portion of the a, x-plane covering this range has been divided into 
various regions. In each region, a distinct method of computation for U(a, x) is used. 
Following a general description, the different regions will be discussed individually. 


General Description of Methods of Computation of Ufa, x) 

For small values of |a| and x, direct computation of U(a, x) by series gives satis- 
factory results. The region containing such combinations of a and x is denoted by S 
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on the figure. The rapidity of convergence of the series decreases as |a| increases, so 
that more terms must be used to calculate a given U(a, x) to a given accuracy. For these 
values of a, the recursion relation for Weber functions can be used to compute U(a, x) 
with the same accuracy as the series computation but in a shorter time. The region 
where this condition prevails is denoted by R g . For certain combinations of |a| and 
x, asymptotic expressions exist that can be used to obtain accurate values for U(a,x) in 
a shorter time than is required by the other procedures. Region Aj denotes the case in 
which 4 1 a | - x^ » 1 and region A n denotes the case in which x^ - 4|a| » 1. Un- 
fortunately, there exists a region for which all the aforementioned techniques, singly or 
in combination, fail to give reliable accuracy. This region is denoted by Rj^ and has 
been omitted from the program. The individual regions will be discussed next. 


Region S - Computation of Ufa, x) by Series 

Region S is defined by the boundaries a[-0. 5, -3], 0 =£ x < 5.5. The even and odd 
solutions (denoted by y^ and y 2 , respectively) to Weber's equation (eq. (4)) can be 
expressed in terms of hyper geome trie functions (ref. 4) as 
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_ _ x 2 \ v /a 1. 1. x 2 

y 1 = expl Ji F i(- + -»-> — 

1 V 4/ 1 x \2 4 2 2 




Then, U(a,x) may be written as 


U(a, x) = YZ yi Yl y 2 ( 7 ) 

2 [(a/2) + (l/4)] r ^a + 1 2 [(a/2) - (1/4)] r ^a + 1 j 

The series expressions for yj and y 2 can be conveniently written in the form 


QO 

'l W * A A 2n« 


where 


A q (x) = 1 


A 9 (x) = a — 
* 2! 


when n ^ 1, and 


l2n+2<X) ■ (2n + lH2n + 2) [ aA 2» W * 4 


QO 

y 2 lx) = E a 2iw1 (x) 


where 
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A i (x) = x 


Aq(x) = a — 
J 3! 


J 


( 12 ) 


and 


"2n+3 


(x) = 


(2n + 3)(2n 


7^ l^Si+lW + 7 x2A 2n-l <x) ] 


(13) 


when n ^ 1 . 


Region R s - Extension of Series Computations by Use of Recursion Relations 

2 

This region is defined by - (50 + x )/4 < a ^ -3.0, 0^x^ 5.5. The recursion 
relation used may be written 

U(a-1, x) = xU(a, x) + ^ + i^U(a+l, x) (14) 

and is called the backward recursion relation. From it, a Weber function of order 
a - 1 is obtained from two Weber functions of higher order, that is, a and a + 1. The 
use of the recursion formula still requires that two values of U(a, x) be computed by 
series, that is, U(a, x) and U(a+l,x). 

It was found that, starting with 8 -figure accuracy, equation (14) can be used 30 suc- 
cessive times with a loss of accuracy only in the last place. If the corresponding "for- 
ward" recursion relation (which can be used to obtain higher orders of U(a,x) starting 
with lower orders) be used, loss of accuracy occurs after the very first use of the re- 
cursion formula. However, using double precision, the forward recursion relation may 
be used about 15 times and still maintain 6 -figure accuracy. 


Region Aj - Asymptotic Solutions for |a| Large and x Moderate 

Region Aj is defined as that region for which 4 1 a | - x^ > 50. In this region a 
function which will be called asymptotic I is used for computing U(a, x). This function 
may be written as 
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where 


U(a,x) ~ 2- 


S3 

y/2^Y 


e r cos/ — + — + 9 + v. 

U 2 1 


f 

2 X 


xY 


Y dx = — + a| sin 
4 


, -1 


L 2 (| a |) 1 /2 


=Vj4| 


a - x 


V 

r y 6 ^12 


v = ^ + 
Vi ^3 " y 9 


j l/x 3 1 J 
d, = —I — + — ax 

d a 1 48 2 > 


d fi = 1 x 2 - 2a 
6 4 


d ft = A/l _1_ x 9 - A_ ax 7 - — a 2 x 5 + A a 3 x 3 - 19a 4 x ' 


.31 5760 320 320 
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d „ = A 3 x 4 - 186ax 2 + 80a 2 
11 8 


(15) 


(16) 


(17) 

(18) 

(19) 

( 20 ) 


This asymptotic approximation truncating v r and Vj after two terms gives 7 -figure 
accuracy if 16 figures are used in the calculations. It should be noted that the form 
of asymptotic I differs from that in reference 4 (p. 690) because of a slightly dif- 
ferent definition of v^. 
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Region Ajj - Asymptotic Solutions for Large x 

At the opposite extreme, where x 2 - 4|a| » 1, a different approximation exists 

which gives satisfactory results. Region A n is defined as that region for which 
2 ii ^ 

x - 4 1 a I > 50. The asymptotic expression used to compute U(a,x) in this region will be 
called asymptotic n and can be written as 

/l/- - a) 

U(a,x)~ 1/— -exp[-0+ v(a,x)] (21) 

where 



r 


Ydx = -xY+aln 
4 


x+ Y 

2(| a |) l/2 


v(a,x)=Z (-l) s — 

S=1 y3s 


and Y, dg, dg, dg, and djg are given by equations (16) to (20). The accuracy of asymp- 
totic n is the same as that of asymptotic I and the same note applies also. 


Region R A ^ - Asymptotic II Plus Recursion 

Region ^Ajj * n figure 1 (p. 4) is a region in which none of the conditions prevail 

which would enable the use of any of the preceding four methods of computation to give 

accurate results. It is essentially the region between Aj and Ajj but above x =v£o7 

In RAtt’ vuiu® 5 °f x are too large to enable series computations to be made ac- 
II 2 i i 

curately. On the other hand, they are not so large that x - 4|a| > 50. Neither is a 

so large that 4|a| - x 2 > 50. Fortunately, the use of the recursion relation (14) applied 

to starting values of U(a, x) in region A n can carry computations through region ^Ajj 

with acceptable accuracy (7 significant figures) . 


Region R. -Asymptotic I Plus Forward Recursion 
A I 

This region is bounded by x = [5. 5, V60], - [(x 2 + 2)/4] > a > - [(x 2 + 50)/4]. 
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Forward recursion from Aj will still give fairly accurate results (5 figures) up to 
15 recursions. Therefore, (a, x) combinations in RAj can be computed by locating the 
point in Aj and. then recurring in the forward direction. The smaller the number of re- 
cursions, the greater the accuracy of the results. 


Region R NG 

This region is bounded on the bottom by x = 5. 5, on the right by a = - [(x + 2)/4], 
by a = -0. 50 between x = 5. 5 and V"5®> by x = ^56" between a= -0.50 and a = -1.50, 
by a = [ (x^ - 50)/4] between x = VEfF" and x = V60, and by x = -y/60 between a = -2. 5 
and a = -15. 5. For x > 5. 5, the accuracy of the series begins to decrease rapidly, and 
by x ® 8, a = -0. 6, only the magnitude of U(a, x) can be relied on from the series com- 
putations. As explained in the preceding section, the accuracy of forward recursion from 
region Aj is also not very great after many recursions have been made. If values of 
U(a, x) for (a, x) combinations in Rj^-q are nevertheless required, the fact that the ac- 
curacy decreases as x increases and as |a| decreases should serve as a guide. 

PROGRAM DESCRIPTION 

The subroutine WEBER computes the Weber function and its derivative for any point 
in the region a[-l. 5, -200. 5], x[0, 30]. The subroutine is referenced by the statement 
CALL WEBER (XD, AD, V, UP, ISI, IFLAG). The first four arguments of the calling 
vector are double precision and the last two are integers. These arguments are defined 
as follows: 

INPUT 

AD a coordinate 

XD x coordinate 

OUTPUT 

V U(a, x) 

UP U'(a,x) 

TOT 

ISI exponent of scale factor 10 of U and U' 

IFLAG indicator showing whether scaling was used or calculation was omitted 
IFLAG = 0, no calculation made 
IFLAG = 1, no scale factor used 
IFLAG > 1, scale factor is 10*®* 

The structure of the subroutine WEBER is based on locating a specified a,x- 
combination in one of the regions of figure 1. The methods of computing U(a, x) in each 
region are those described in the previous section. To compute U’(a,x) both U(a,x) and 

9 



U(a+l,x) must be available (eq. (3a)). In regions S, Aj, or A^, U(a+l,x) must be com- 
puted separately, but in regions Ra^ r Ajj> and R g> the value comes as a byproduct of 
recursion. 

An outline of the sectional breakup of WEBER follows. (Refer to the flow charts and 
listings in appendixes A and B, respectively, for a more detailed description. ) 


block 1 - Constants Defined 

Constants used in WEBER as well as in other subroutines through COMMON are 
computed and given literal names. 


block 2 - Region of (a,x) Determined 

The process by which any input point (a,x) is located in a region of figure l(p. 4) is 
as follows. The case when 5. 5 < x < V£o" is checked to see if the point (a,x) is in the 
special regions Raj or Rjjq- In the event that (a,x) is in Rj^q, that is, 

(4|a| - x 2 < 2, IFLAG is set to 0 and no calculation is made. When 
2 < (4 1 a J - x 2 ) < 50, the point lies in RAj and two base points in region Aj are picked 
at which U is evaluated for use in forward recursion to obtain U(a,x). Program flow is 
transferred to BLOCK 6. For points in the other regions, there is the following locating 
process. If |a| ^ [(x 2 + 50)/4], (a,x) is in region Aj. If |a| ^ [(x 2 - 50)/4], (a, x) is 
in region Ajj. The procedure for the remaining possibility where 
[(x 2 + 50)/4] < | a | < [(x 2 - 50)/4] is as follows: if x s - -^60^ then (a, x) is in region 
R a ; if x^ 5.5, (a, x) is in S if |a| ^ 3, otherwise it is in Rg. 

block 3 - Series Evaluation in Region S 

The subroutine SERIES is called in BLOCK 3 for two adjacent points in region S (the 
same x coordinate but one unit apart in a). The two values of U are stored for use in 
the recursion relation or in the evaluation of U' . 


block 4 - Evaluation by Asymptotic II in Region Ajj 

The Weber function is evaluated for two adjacent points in region Ajj. As in BLOCK 
3 these values are stored for use in either the recursion relation or for evaluating the 
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derivative. The logic within the DO-LOOP chooses the correct gamma function sub- 
routine for the magnitude of a and keeps the same scale factor for both U values when 
a < -50.0. 

For example, if a = -45.0 and x = 17.6, find U(a,x)and U'(a,x). Since 
R = 129.76 and x 2 > |a|, the program locates (a, x) in Ajj. Then 

a^ = a + 1.0 = -44.0 
= aj — 1. 0 = —45. 0 

The DO -LOOP in BLOCK 4 computes U for these two values of a by using asymp- 
totic H. The results are 


U(-44. 0, 17. 6) = 1. 0484715X10 19 
U(-45. 0, 17. 6) = 1. 5357027X10 20 

Coming out of the loop, a is set equal to a 2 , which is the original a coordinate. 
This condition signals the program to use the two computed values to evaluate U'(a, x). 
The result is 

U’(-45. 0,17.6)= -i(17.6) U(-45. 0, 17. 6) - (-44.5) U(-44. 0, 17. 6) = -8. 848485 lXlO 20 
2 

block 5 - Evaluation in Regions and R$ 

Values of UfegjXj) and Ufeg+ljXj) coming from either BLOCK 3 or BLOCK 4 are 
used in the recursion equation to generate the next value U(a 2 -l,x 1 ). When moving 
along the line x = Xj in the negative direction in a, U(a, x) is evaluated at every point 
UCag-mxj), n = l,2 . . . , until U(a^,Xj) is found. Then U' (a, x) is computed. 

For example, if a ^ -40.00, x = 12.00, and R = |(12.0) 2 - 4 1 —40 . 0 1 J = 16.0, find 
U(a,x) and U*(a,x) when (a,x) lies in R^. 

It is noted that the point (-40.0, 12.0) is closer to region Aj than to region Ajj on 
the line x=12.0. Nevertheless, A n base points will be used because backward recur- 
sion is the accurate one. The two base points for recursion are (-22. 0, 12. 0) and (-23. 0, 
12 . 0 ). 

From BLOCK 4 the following results are obtained: 

U(-22. 0, 12.0) = 6. 041521 1X10 6 
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U(-23. 0,12.0) = 5.9467925X10 7 
Applying equation (14) repeatedly gives 

U(-39. 0, 12. 0) = 3. 5549969X10 22 
U(-40. 0, 12. 0) = 2. 2109229X10 23 


Using equation (3a) gives 


U’(-40.0, 12.0) = -6.0 U(-39. 0,12.0) - (-39.5) U(-39. 0, 12. 0) 


= 7. 7670203X10 22 

block 6 - Evaluation by Asymptotic I in Region Aj 

In this DO -LOOP as in BLOCK 4, U is evaluated at two adjacent points. These 
points (a,x) and (a+l,x) are used in computing U' . The same type of coding as in 
BLOCK 4 is used for finding the gamma function and keeping the scale factor consistent. 
In case the calculations are made to provide base points for recursion, the following 
forward recursion relation is used until U is evaluated for the desired a: 

U(a x) = U(*-2,x) - xU(a-l,x ) 

a. - — 

2 


The derivative is evaluated by the relation 

U’(a,x) =ixU(a,x) - U(a-l,x) 

2 

Subroutines 

Of the seven subroutines, five are used for evaluating the gamma function. Sub- 
routines GAMMA, GAMMAF, GAMMAN, STRGAM, and STRSCL have been developed to 
evaluate the gamma function in the range - 2 < z < 200. The two approximations used 
are (1) the series expansion to 26 terms for l/r(z) and (2) Stirling’s approximation 
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r(z) ~ z ( z -l/2) e - z ( 2irP^(l + — + — 

y 12z 


139 


571 


288z 2 51840z 3 2488320z 4 


+ . . . 1 ( 22 ) 


Subroutine GAMMA is written to evaluate r(z) by computing 


26 


Cyf for 1 < z < 2 


and then taking the reciprocal. The C k are displayed in the listing for GAMMA. The 
maximum error of r(z) in this range of argument is 2xl0~ 3 when 16-place numbers are 
used. GAMMA is called by two subroutines; either GAMMAF, which monitors the argu- 
ment z s 1 , or GAMMAN, which handles nonintegral arguments z < 1. GAMMAF re- 
duces the positive argument to within the range of GAMMA and then uses the relation 


T(z + 1) - zr(z) (23) 

GAMMAN, on the other hand, increases the argument until it falls within the range used 
in GAMMA and uses the relation (eq. (23)) in the negative direction; that is, 

r(z) = r k t Jl (24) 

z 


STRGAM and STRSCL are altered forms of Stirling’s approximation 
r(z), where z = 1/2 - a, appears in both asymptotic I and n in the form 
equation (20) can be rewritten to compute this factor directly by using 


(e q. (2 0)). Since 
|/r^)/(2 7 r) 1 / 4 , 


Vr(z) a /z\ z / 2 z -l/4 L + J_ + 1 _ 139 

(27T) 1 / 4 \ 12z 288z 2 51840z 3 


571 

2488320z 



1/2 


(25) 


STRGAM computes this factor for | a | > 10. 0 with at least 8-figure accuracy when 
16 figures are used. 

STRSCL computes the same factor for |a| >50.0 with the added feature of a scale 
factor, since the value of the altered gamma function becomes too large for the computer. 
Scaling is done on equation (25) as follows: Let 


13 



and 


y = 



KI = [y] 

greatest integer function. Then ESI, a power of 10, becomes the scale factor and 
lO.ofr- 181 ) replaces (z/e) Z/ ^ in equation (25). The ISI is carried back to WEBER in 
the calling vector STRSCL. 

Scaling of the gamma function is meshed with scaling of the factor exp[-0 + v(a., x)] 
in asymptotic n. Whenever the argument [-0 + y(a, x)] approaches the value -88, below 
which the exponential would be set to zero, a scale factor is factored from the exponen- 
tial and combined with the scale factor of gamma. The quantity 

exp £-0 + i/(a, x)J is written 10 {[ ^ a,x ^] ^lO 6 } 


The exponent of 10 is split into the largest integer and a remainder. The integer is com- 
bined with S coming from STRGAM to give a scale factor for U, and the value of 10 
raised to the remainder replaces exp [-6 + v(a., x)] in the computation of U. 

The other two subroutines are SERIES, which computes U(a, x), equations (7) to (13) 
andDFIND, which evaluates the coefficients dg^, equations (17) to (20). 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, October 5, 1966, 

125-23-02-10-22. 
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APPENDIX A 


FLOW CHARTS FOR COMPUTATION OF WEBER FUNCTIONS AND FIRST DERIVATIVES 
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APPENDIX B 


PROGRAM LISTINGS FOR COMPUTATION OF WEBER 


FUNCTIONS AND FIRST DERIVATIVES 

MAIN DEBUG 

DOUBLEPRECISION A, X. U. UP* S , AMI N. AMAX.DA. XMAX. DX. S AVE. DABSU. DABSUP 
AMIN * -5. DO 
AM AX * —200. DO 
DA * -5. DO 
XMAX * 30. DO 

- DX = I. DO 

A = AMIN 

\ 25 WRITE 16.1) A. A 

1 FORMAT (1H1 6X 1HX 7X 2HU1 F6. L « 3H.X) 10X 2HU1 F6.1.3H.XI/1HJ) 

X * O.DO 

ILJA = 0 

50 CALL WEBER I X , A. U. UP . I SI . I SCALE ) 

IF 1 ISCALE.EO.O) GO TO 120 
IS = 0 
JS = 0 

OABSU * DABS ( U ) 

IF 1 DABSU ) 60.175.60 
60 IF ( DABSU.LT. 1 .DO I GO TO 300 
IF 1DABSU.LT. 10. DO) GO TO 100 
C NUMBER GT 10.0 
200 U = U* 1 *D— 1 
IS = IS + 1 
DABSU = DABSIU) 

IF ( DABSU -10.DOI 100.200,200 
C NUMBER LT 1.0 
300 U = U* 1.D1 
IS = I S-l 
DABSU = DABSIU) 

IF 1 DABSU— 1.00) 300.100,100 
100 DABSUP * DAB SI UP ) 

IF 1 DABSUP.LT • 1. DO ) GO TO 400 
IF IOABSUP.LT. 10. DO) GO TO 700 
C NUMBER GT 10. 0 
500 UP = UP* 1 .D— 1 
JS = JS+1 
DABSUP = DABS! UP I 
IF (DABSUP-iO.DO) 700.500,500 
400 UP = UP* 1 • D+l 
JS - JS-1 
DABSUP * DAB SI UP ) 

IF ( DABSUP— 1 .DO ) 400,700.700 
700 IF IISCALE.EO.il GO TO 75 
IS = IS+ISI 
JS = JS+ISI 
75 XX = X 

uu = u 

UUP = UP 

* IF (ILJA.EO.OI GO TO 118 

IS = IS- I SI 
JS = JS-JSl 

WRITE <6.21 XX.UU.IS.UUP.JS 

2 FORMAT <5X F4.1, 3X F10.7. IX 1H( I4.1H) 3X F10.7.1X IHC I4.1H)) 

; GO TO 120 

118 WRITE (6,3) IS,JS 

3 FORMAT 1 24X 14. 16X 141 

JS1 = JS 

IS1 * IS 
ILJA = 1 
JS s 0 
IS = 0 

WRITE (6.2) XX. UU. IS, UUP, JS 
120 X = X+OX 

IF ( X— XMAX ) 50.50.125 
125 IF IA-AMAX) 250,250.150 
150 A = A+DA 
GO TO 25 
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175 XX * X 

uu = u 

UUP - UP 

WRITE 16*2) XX.UU*IS,UUP.JS 
GO TO 125 
250 STOP 
END 


SUBROUTINE WEBER ( XD*AO* V* UP* I SI * I FLAG) 

C FOR THE REGION A I- 1 .5,-200. 0) * X 10-0*30.0) 

C THE CALLING VECTOR 
C XD - THE X COORDINATE 
C AD - THE A COORDINATE 
C V - U( A* X I 

C UP — DERIVATIVE OF U(A*X1 
C ISI - SCALE FACTOR 

C I FLAG = 2 OR 3 WHEN 10.0**S IF USED AS A SCALE FACTOR OF V ANO UP 
C I FLAG = 0 WHEN NO CALCULATION MADE 

COMMON CA3* CB2,CC9,CC7.CC5*CC3*CCl .C04.C02.C00 
DOUBLE PRECISION A* ASA VE . AB SA. ARG* ANS . ANGLE .CHI ,CHI 0 

* *CA3*CB2 « CC9* CC7.CC 5*CC 3. CC1 .CD4.CD2 «CDO .03* 06 *D9 .012 

* * GAM* PI*PI2*PI 24 »R * S* THE TA* U* UP • V* VAX * VI * VRP. X* X2 *XD* AO 

* * SAM* TAM* DLOGE*RP 
DIMENSION I S( 2 ) * Ut 3) 

I FLAG = 1 

CCC BLOCK 1 - SETUP DATA 
X - XD 
A = AD 

PI = 3.141592653589793 
P 12 = 6. 283185307179586 
P 1 24 = PI2**.25 
C COMMON CONSTANTS 

CA3 * 1.D0/24.D0 
CB2 = 75. D- 2 
CC9 = 7. D0/5760. DO 
CC7 = 7. 00/320. DO 
CC5 = 49. D0/320. DO 
CC3 = 31. DO/12. DO 
CCl = 19. DO 
CD4 = 153. DO/8. DO 
CD2 = 186. DO 
CDO = 80. DO 

DLOGE - .4342944819032518300 
ASAVE = A 
X2 = X*X 

CCC BLOCK 2 - DETERMINE REGION OF (A. XI AND EVALUATION METHOD 
ABSA = DABS(A) 

IF 1X2 —60. DO ) 2.4*4 

2 IF 1X2 -56. DO) 3*5.5 

3 IF (X-5.5D0) 6*6.7 

4 IF (ABSA . GE • ( 5 0.D0+X2) /4. DO) GO TO 200 

11 IF (ABSA *LE.(X2— 50. DO) /4. DO) GO TO 30 

A = A +1.D0 

ABSA = DABS! A I 
GO TO 11 

5 IFt ABSA.LE.l X2— 50.D0 ) /4.D0 ) GO TO 30 
13 IF( ABSA.GE. ( 50.00 + X2)/4.D0) GO TO 200 

IF (ABSA .GE. ( X 2+2. DO) /4. 00) GO TO 12 
I FLAG * 0 
WRITE! 6. 161 

16 FORMAT ( 1H0 2X 23HP0INT IS IN BAO REGION / ) 

RETURN 

12 A = A-l.DO 
ABSA = DABS(A) 

GO TO 13 

6 IF ( ABSA.GE. ( 50. D0+X2I/4.D0) GO TO 200 
9 IFt ABSA.LE. 3. DO) GO TO 20 

A * A+l.DO 
ABSA = DABS t A ) 

GO TO 9 

7 IF( ABSA . G€ • ( 50.D0+X2) /4.D0I GO TO 200 
IF (ABSA • GE. ( 2 .D0+X2 ) /4. DO) GO TO 8 

I FLAG * 0 
WRITE (6*17) 

17 FORMAT ( 1H0 2X 23HP0INT IS IN BAD REGION / ) 

RETURN 

8 A = A-l.DO 
ABSA * DABS(A) 

GO TO 7 
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CCC BLOCK 3 - SERIES EVALUATION IN REGION S 
20 A - A+I.DO 
DO 25 N*l,2 
CALL SERIES (X,A,ANSI 
UC M I * ANS 
25 A * A— 1.00 
A * A+l.DO 

IF I A— ASA VE I 60,60,50 
CCC BLOCK 4 - ASYMPTOTIC II IN REGION All 

30 IF ( ASAVE.LT. -50. 00) IFLAG * 2 
ISCl) * 0 

ISC 2) * 0 
A * A+l • DO 

00 40 M*l, 2 
ABSA * OABSCA) 

CHI * DSORTC X**2-4.D0*ABSAI 

CHID * l.00/CHI**3 

CALL OFINOI X, A, 03, 06,09,0 12) 

VAX * -,5DO*DL0GCCHI)+CHID*(-O3+CHID*(D6+CHID*(-D9+CHID*DI.2) ) 1 
THETA » .2500* X*CHI-*-A*DLQG( ( X+CHI ) /I 2* D0*DS0RT C ABSA) ) ) 

ARG » • 5 DO- A 

IF I ABSA— 10*001 31, 31,32 

31 CALL GAMMAF I ARG, GAM I 

C POLYNOMIAL APPROXIMATION OF THE GAMMA FUNCTION 
GAM * DSQRKGAMI/PI 24 
U< Ml * GAM*DEXPC-THETA+VAX> 

GO TO 40 

32 GO TO (36,33,36), IFLAG 

33 CALL STRSCL (ARG, GAM, SI 

C SCALED, REOUCED STERLINGS APPROXIMATION OF GAMMA FUNCTION 
IS(M) = S 

SAM * C-THETA+VAXl*DLOGE 
IL * I DINT I SAM ) 

TAM * SAM— OBLE ( FLOAT ( IL I > 

IS(M) = I S( Ml +IL 
IF I M— 2 I 34,35,34 

34 U(NI * GAM* 10.D0**TAM 
GO TO 40 

35 N = IS( 2 1 - ISI 11 

TAM = OBLE(FLOAT(N) H-TAM 
U(MI * GAM* 10 -DO** TAM 
ISI * ISC II 
GO TO 40 

36 CALL STRGAM (ARG, GAM) 

IF l IFLAG, EQ.3.AND.M.EQ.2) GO TO 362 
SAM = -THETA+VAX 

IF (M.EO. I.AND.SAM.LT.-88.D0I GO TO 361 
U(M) = GAM *0E XP ( SAM ) 

GO TO 40 

361 IFLAG * 3 

SAM = C— THETA+VAX )*DLQGE 

IL = I D INK SAM I 

TAM = SAM— DBLE (FLOAT! IL) I 

UCM)=* I0.D0**TAM* GAM 

ISC M ) = IL 

GO TO 40 

362 SAM * (—THE TA*VAXI*DLOGE 
IL = I D INK SAM ) 

TAM * SAM— DBLE ( F LOAT ( I L ) ) 

1 SC M ) * IL 

IM = IS( 21— IS( I > 

U( M I - GAM*10.D0**TAM*10.D0**IM 
ISI * ISC 1) 

40 A * A— 1,00 
A * A+l.DO 

IF (A-ASAVE) 60,60,50 

CCC BLOCK 5 - RECURSION EQUATION FOR REGIONS RA AND RS 
50 UC 3) * (A+.5D0i*UC 1) + X*UC2) 

U( I) * UC 2) 

UC 2 ) * UC 3) 

A * A- 1,00 
DEBUG A, U( 2 ) 

IF C A-ASAVE) 60,60,50 
60 UP * — • 5D0*X*U( 2 I - (A+.5D0)*U41I 
V * UC 2 1 
RETURN 



CCC BLOCK 6 - ASYMPTOTIC I IN REGION AI 
200 IF I A.LT.-5O.D0) IFLAG * 2 
A * A+l-DO 
00 240 M*l* 2 
ABSA s OABS(A) 

CHI * OSORTI 4*D0*ABSA-X*XI 

CHIO * l.D0/CHl**3 

CALL OF I NO I X* A, D3, D6*D9,D12 ) 

VRP * CHID*CHID*(-D6*D12*CHtD*CHIDI 
VI = CHID*! 03— D9*CH ID*CHIO ) 

THETA = *25D0*X*CHI+ABSA*DATAN IX/CHI ) 

ANGLE * ( THETA+VI *PI *1 -25D0+* 500*A ) I 
ANGLE * DMODI ANGLE * P I 21 
ARG * -500- A 
GO TO I 205*210) * IFLAG 
205 CALL STRGAHI ARG* GAM ) 

UI M ) * 2-D0*GAM*DEXPI VRP) /DSQRTICHI)*DC0S1 ANGLE) 

GO TO 240 

210 CALL STRSCL ( ARG*GAM* S) 

IS! M ) * S 

IF ( M— 2 ) 220*230,220 

220 UIH) = 2 *DO*GAH*DE XP( VRP ) /OSORTI CHI )*DCOSIANGLE) 

GO TO 240 

230 N * ISI 2 )— I SC 1) 

UIM) * 2-00*GAM*10-00**N*DEXP( VRPI /OSORTI CHI )*DCOS I ANGLE ) 

is! = ism 

240 A = A— 1 • DO 
A = A+1-00 

IF CA-EO-ASAVE) GO TO 260 
C FOREWARO RECURSION FOR REGION RAI 
A * A+l.DO 

245 UI3) « IUI2I - X4UI1II/IA+.5D0) 

A = A+l-DO 

IF ( A-EO-ASAVE ) GO TO 250 
UI2J = UI 1) 

Vi 1) * UC3I 
GO TO 245 

250 UP = • 5D0*X*UI 3 ) - UCII 

V = UI 3) 

RETURN 

260 UP * — • 5D0*X*UI 21— I A+-5D0 ) *UI 1 ) 

V = UI 2 I 
RETURN 
END 


SUBROUTINE SERIES IX*A,ANS) 

C COMPUTES UIA.XI BY EVALUATING SERIES Yl AND Y2 
DIMENSION CZEROm,Cf 10001 

DOUBLE PRECISION A, ARG* ANS*C,CZERO*C YQNE *C YTWO * ERROR *FN* GAM 
* * PI 2, TW0A«TW04*X, X2* YONE * YTMQ*CARG* SONE*STWO 
PI2 = 1.772453850905516 
TW04 * 1*1 892071 l 5D0 
ERROR = 1*D-10 
SONE = O-DO 
STWO * O.DO 
X2 = X*X 
CZERO * 1.00 
Cl 1) * X 

C ( 2 ) * A4X2/2-D0 
Cl 31 = A*X2*X/6*D0 
YONE * CZERO *01 21 
YTWO = CC1I+CI3I 
DO 100 NT * 2*900,2 
FN * DBLEI FLOAT! NT) ) 

CtNT+2) « X2/CIFN + l*D0)*<FN+2. DO) 1 *< A*CI NT ) ♦• 25D0*X2*C C NT— 2 I ) 
CINT+3) * X2/I I FN+2-D0)*! FN+3.D0) ) *1 A*CC NT* II ♦•25DQ*X2*C( NT— 1 ) ) 
YONE * YONE+CI NT+2 ) 

YTWO * YTWO+CCNT+3) 

IF I DABS! YONE- SONE) - ERROR) 99,99,101 
99 IF I DABS! YTWO- STWO) - ERROR) 110,110,101 



101 SONE * YONE 
STWO * YTWO 
100 CONTINUE 
110 ARG * • 5DO*A 

TWOA * 2.00**ARG 
ARG * .75D0+.5D0*A 
CARG * OBL El FLOAT! I DINTt ARG) 1 1 
IF I ARG-CARG I 130,120,130 
120 CYONE * O.DO 
GO TO 140 

130 CALL GAMMANI ARG,GAM) 

CYONE * P12*Y0NE/ITW04*TW0A*GAMI 
140 ARG = . 25D0+. 5D0*A 

CARG = DBLE (FLOAT! IOI NT! ARGI I ) 

IF ( ARG— CARG I 160,150,160 
150 CYTWO * 0.00 
GO TO 170 

160 CALL GAMMAN I ARG, GAM I 

CYTWO * PI26TW04*YTW0/(TH0A*GAMI 
170 ANS * C YONE— C YTWO 
RETURN 
END 


SUBROUTINE DF IND( X, A , D3 , D6 » 09, D1 2 ) 

COMMON CA3, CB2 , CC9 , CC7, CC 5 , CC3 , CC 1 , CD4 , C02 , COO 

DOUBLE PRECISION X , A , D3 , D6 , D9 , D I 2 , A2 , X2 , C A3 , CB2 , CC9 , CC7 , CC5 , CC3 

• , CCI,CD4,CD2» CDO 
A2 = A*A 

X2 = X*X 

D3 * .500*X*(CA3*X2+A)/A 
D6 » CB2*X2-2 • DO*A 

09 = < ( ( (-CC9*X2-CC7*A)*X2-CC5*A2)*X2 + CC3*A2*A)*X2- CC1 *A2*A2)*X 

• / ( A2*A) 

012 = X2* ( CD4* X2-CD2 *A )+CD0*A2 

RETURN 

END 


DOUBLE PRECISION FUNCTION GAMMA ( X ) 

DOUBLE PRECISION Z , C 1 , C 2 , C3 , C4 , C5 , C6 , C 7 , C8 , C9 , CIO , C l 1 , C 1 2 , C 1 3 , C 14, 
A C15#C16fC17,C18,C19 t C20,C2l,C22,C23,C24,C23,C26,RECIPR,X 
DATA C1»C2,C3»C4,C5, C6,C7,C8,C9,C10,C11 ,C12,C13,C14,C15,C16,C17, 

A C18»C19,C2Q»C21,C22»C23,C24»C25»C26 / 1 . OCOOOOOOOOOOOOOODO , 

A 0. 5772 1 566490 1 5329D0 , -0. 655878071 5202 5 38D0, -0.042002635034095200, 
B 0. 166538611382291500,-0.042197734555544300,-0.009621971527877000, 
C 0. 0072 1 894324666 30D0, -0.001 165 167591 8 59 100, -.00021 5241 674 114900, 

0 0.000 12 8050282 388200,-0. 000020 134854780 700, -0.00000 12 50 49 3482 IDO, 
E 0.00000 11 33027232000,-0. 000000205633841 7D0 , 0. 0000000361 16095000, 

F 0. 00000000 500200 7500, -0.00000000 l 18 12746D0, 0.000000000 104342700, 

G O.OOOOOOOOOOOT7823DO, -0. 000000000003696 8D0, 0. OOOOOODOOOOD510DDO, 

H -0.000000000000020600,-0. 000000000000005400, 0.00000000000000 1400, 

1 O.OOOOOOOOOOOOOOOIDO/ 

Z * X 

IF (Z.LT.1.DO.OR.Z.GT.2.DO) STOP 

RECIPR = Z*(Cl+Z*(C2+Z*(C3+Z*(C4+Z*(C5+Z*(C6+Z*(C7tZ*(C8+Z*(C9+ 

A Z*(C10+Z*(C11+Z*(C12+Z*(C13+Z*(C14+Z*(C15+Z*< C16 + Z*IC17* 

B Z*(C18+Z*<C19+Z*(C20+Z*(C21+Z*(C22+Z*(C23+Z*(C24+Z*(C25+ 

c z*c 26 1 m n m m m m m m ) 

GAMMA = l.ODO/RECIPR 

RETURN 

END 
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SUBROUTINE GAMMAF ( X , GAM ) 

C FOR POSITIVE VALUES OF X GREATER OR EQUAL TO 1.0 
DOUBLE PRECISION . X, Z , ARG, GAM , GAMMA 
Z = X 
ARG = Z 

10 IF (ARG-2.D0) 20,20,15 
15 ARG = ARG- 1 • DO 
GO TO 10 

20 GAM = GAMMA (ARG) 

30 IF (ARG-Z) 35,40,35 
35 GAM = GAM*ARG 
ARG = ARG+l.DO 
GO TO 30 
40 RETURN 
END 


SUBROUTINE GAMMAN ( X * GAM ) 

C FOR NON-INTEGRAL VALUES OF Z RANGE (-33.0,1.0) 
DOUBLE PRECISION X, Z , ARG, GAM, COEF, GAMMA 
Z = X 

IF (Z) 10,5,5 
5 ARG = Z+1.D0 

GAM = GAMMA ( ARG ) /Z 
RETURN 

10 COEF = i .00 
20 COEF = Z*C0EF 
Z = Z+1.D0 

IF (Z-1.D0) 20,20,30 
30 GAM = GAMMA ( Z ) /COEF 
RETURN 
END 


SUBROUTINE STRGAM (Z,ANS) 

C COMPUTES SQRT ( GAMMA ( Z ) ) / ( 2. »P I ) ** . 25 

DOUBLE PRECISION Z , ANS , A, B , C , D, E , SER I ES , Z1 , POW 
DATA A,B,C,D/83333333333333333.D-1 8, 34722222222222222. D-19, 
A 268 1327 1604938272. D-19, 22947209362139918. D-20/ 

DATA E 72.7182818284590452/ 

POW = Z/2.OD0 
Z1 = l.ODO/Z 

SERIES = 1.0D0+Z1*(A+Z1*(3-Z1*(C+ZZ*D) ) ) 

80 ANS = (Z/E)**POW/Z**.25DO*DSQRT( SERIES) 

RETURN 

END 


SUBROUTINE STRSCL (Z,R,S) 

C COMPUTES SQRT ( GAMMA (Z))/(2.*PI)**.25 WITH SCALING FACTOR S AS 
C A POWER OF 10.0 

DOUBLE PRECISION Z, Z1 , A, B, C, D, E, SERIES, S , P, ZE, R 
DATA A, B,C, D/83333333333333333. D-18, 34722222222222222. D-19, 
* 268 13271604938272. D-19, 2294720936213991 8. D-20/ 

DATA E /2. 7182818284590452/ 

Z1 = l.ODO/Z 

SERIES = 1.0D0+Z1*(A+Z1*(B-Z1*(C+Z1*D) ) ) 

S = C.00 

ZE = DL0G10 ( Z/E ) 

P = Z/2.D0*ZE 
10 S = S+l.DO 
P = P-1. DO 

IF (P-1. DO) 20,10,10 
20 R = 1O.D0**P*DSURTISERIES>/Z**.25DO 
RETURN 
END 
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APPENDIX C 


TABLES OF THE WEBER FUNCTION AND ITS FIRST DERIVATIVE 


The following tables list D^(x) (U(a, x) column heading) and (d/dx)D^(x) (U'(a, x) col- 
umn heading) at intervals of 1 in x and 5 in v . The true exponent is obtained by adding 
the number at the top of each exponent column to the exponent of each individual entry. 


X 

U< -5.0 

• XI 

UM -5. 

O.X) 



0 


0 

0. 

3.0521837 

( Oi 

6.8415 762 

* 0) 

1.0 

5.7992601 

I -1) 

-9.4555887 

< 0) 

2.0 

-4. 1865842 

< 0) 

3.3663954 

i 0) 

3.0 

3.2021291 

( 0) 

6.8170127 

( 0) 

4.0 

5.0331933 

( 0) 

-2.3677349 

I 0) 

5.0 

1.8799768 

< 0) 

-2.6935859 

( 0) 

PO INT 

IS IN BAD REGION 



PO INT 

IS IN BAO REGION 



8.0 

1.1456761 

< -3) 

-3. 8993322 

{ -3) 

9.0 

2. 8559484 

1 -5 ) 

-1. 1357113 

{ -4) 

10. 0 

4.0491466 

( -7) 

-1. 8355745 

1 -6) 

11. 0 

3. 3102931 

< -9) 

-1.6811175 

( -8) 

12.0 

1.5756262 

( -11) 

-8. 8479091 

( -11) 

13.0 

4.3974120 

( -14) 

-2. 7028308 

( -13) 

14.0 

7.2347562 

C -17) 

-4. 82 74995 

I -16) 

15.0 

7. 0459403 

( -20) 

-5. 06 96986 

t -19) 

16.0 

4.0754359 

< -23) 

-3. 1441223 

< -22) 

17.0 

L. 4037342 

I -26) 

-1.1555568 

( -25) 

IB.O 

2.8854817 

i -30) 

-2.5240025 

( -29) 

19.0 

3.5461783 

( -34) 

-3.2840529 

( -33) 

20.0 

2.6095865 

i -38) 

-2.6095865 

I -37) 

21. 0 

0 . 

t 0) 

-0. 

< 0) 






X 

Ul -15.0.X) 

UM -15. 

0 , X) 



5 


5 

0 . 

-1. 8569056 

( 0) 

-7.1937571 

( 0) 

1.0 

2.6317424 

( 0) 

6.6710645 

t -1) 

2.0 

-2.2158064 

( 0) 

5.5473151 

C 0) 

3.0 

1.202761 8 

< o) 

-8.7388111 

< 0) 

4.0 

-4.0493298 

i -l) 

9.3030182 

( 0) 

5.0 

3.9777795 

t -i) 

-8. 7929832 

( 0) 

6.0 

-1.7039120 

< 0) 

6. 7417567 

i 0) 

7.0 

3.9126172 

< 0) 

8. 8378327 

( -l) 

8.0 

1.8565878 

I 0) 

-2.740131 2 

( 0) 

9.0 

2. 5128416 

i -1) 

-6.2181926 

< -1) 

10.0 

1.4010081 

( -2) 

-4.5933533 

( -2) 

11.0 

3.6757464 

1 -4) 

— 1 . 46 71 97 0 

( -3) 

12.0 

4.8632681 

t -6) 

-2.2624029 

( -5) 

13.0 

3.3847003 

1 -8) 

-l. 7866440 

( -7) 

14.0 

1. 2746504 

( -10) 

-7.4970285 

I -10) 

15.0 

2.6501545 

( -13) 

-1. 7139929 

C -12) 

16.0 

3.0877971 

( -16) 

-2.1739371 

( -15) 

17.0 

2.0394185 

i -19) 

-1.5506080 

I -18) 

18.0 

7. 7 C49507 

{ -23) 

-6.2856161 

( -22) 

19.0 

1.6772477 

I -26) 

-1.4602201 

l -25) 

20.0 

2.1162772 

( -30) 

-1. 9573032 

< -29) 

21.0 

1.5553965 

< -34) 

-1.5222734 

< -33) 

22.0 

6.6866087 

1 -39) 

-6.901562 2 

l -38) 

23.0 

1.6873243 

{ -43) 

-1.8311915 

< -42) 

24.0 

2.5068728 

( -48) 

-2.8530711 

( -47) 

25.0 

2.1985814 

i -53) 

-2.6178340 

< -52) 

26.0 

1. 1408229 

1 -58) 

-1.4181251 

i -57) 

27.0 

3.5093471 

( -64) 

-4.5455334 

< -63) 

28.0 

6.4110844 

« -70) 

-8.6376010 

( -69) 

29.0 

6.9664549 

I -76) 

-9.7472657 

( -75) 

30.0 

4.5089107 

C -82) 

-6.5420684 

1-81) 


X 

c 

1 

0 

1 

o 

♦ X) 

U’ l -10. 

0 

• X) 



2 



3 

0 . 

-3. 7799709 

( 0) 

1.1960745 


0) 

1.0 

3. 7675627 

I 0) 

-1.1949964 


0) 

2.0 

-4.1101075 

( 0) 

1.0804553 


0) 

3.0 

5.0326142 

( 0) 

-7. 1888978 


-l) 

4.0 

-6.0616342 

( 0) 

-7. 8046901 


-2) 

5.0 

4.0529365 

I 0) 

1.1308804 


0) 

6.0 

6.5793091 

( 0) 

-4.3712358 


-1) 

PO INT 

IS IN BAD REGION 




8.0 

2.0577253 

t -1) 

-5.3459455 


-2) 

9.0 

1.0686247 

( -2) 

-3.5316828 


-3) 

10.0 

2.8258753 

I -4) 

-1.1171334 


-4) 

11.0 

3.9826571 

1 -6) 

-1.8185412 


-6) 

12.0 

3.0790628 

( -8) 

-1.5874445 


-8) 

13.0 

1.3318242 

( -10) 

-7.6294212 


-11) 

14.0 

3.2690072 

i -13) 

-2. 0559890 


-13) 

15.0 

4.6018429 

( -16) 

-3. 1480659 


-16) 

16.0 

3. 7458023 

( -19) 

-2. 7663513 


-19) 

17.0 

1.7744241 

I -22) 

-1.4060121 


-22) 

18.0 

4.9172834 

( -26) 

-4. 1588659 


-26) 

19.0 

8.0055994 

i -30) 

-7.1951735 


-30) 

20.0 

7.6841260 

I -34) 

-7.3110480 


-34) 

21.0 

4.3612678 

t -38) 

-4.5793312 


-38) 

22.0 

0. 

( 0) 

-0. 


0) 


X 

UI -20.i 

0, 

X) 

U’ ( -20. 

O.X) 




8 


6 

0. 

2. 1958943 


0) 

-9.8218701 

( 0) 

1.0 

1.5923444 


0) 

1.1911176 

( 1 ) 

2.0 

-3.0617271 


0) 

3. 1032601 

< 0) 

3.0 

6.0073639 


-1) 

-1. 3230133 

I 1) 

4.0 

2. 3168492 


0) 

9.3812065 

I 0) 

5.0 

-3.4021392 


0) 

-9.8485444 
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